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In  this  analysis,  the  effect  of  heating  as  a  means  of  decontamination  of  substrates  Is 
examined  in  detail.  In  particular,  the  use  of  thin  electrically  conducting  films  Imbedded  in 
substrates  is  considered  as  a  heat  source.  Solid  substrates,  when  exposed  to  Mndeslrable 
vapors,  can  experience  temporary  or  permanent  material  damage  by  the  adsorption  of  these 
vapors.  However,  the  supply  of  thermal,  energy  to  the  contaminant  can  put  them  back  into  the 
gaseous  state. 

In  the  current  development,  a  one-dimensional  model  for  infinitely  large  substrates  is  , 
adopted.  The  thermal  transport  mechanisms  include  heat  generation  in  the  conductive  layer, 
unsteady  heat  conduction  through  the  solid  substrate,  heat  of  desorption  in  the  adsorbed 
layer,  and  thermal  convection  In  the  gaseous  exterior.  The  mass  transport  of  the  contaminant 
takes  place  by  diffusion  in  the  substrate  and  desorption  at  the  surface.  These  phenomena  are 
mathematical ly  modeled  by  a  set  of  governing  differential  equations  and  boundary  conditions. 
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This  set  of  equations  are  then  made  nondlmensional  by  suitably  chos  ;r.  dimensionless 
groups  and  programed  for  solution  by  finite-difference  methods.  A  listing  of  the 
program  together  with  the  required  Input  parameters’ is  given. 

The  program  is  capable  of  predicting  the  temperature  and  the  contaminant  concen¬ 
tration  of  a  solid  undergoing  heating.  Some  special  cases  have  been  run  to  exhibit 
the  typical  behavior  of  such  systems.  For  most  cases  the  heating  decontaminates  the 
solid.  However,  for  seme  cases  the  heating  increases  the  solubility  of  the  contaminant 
and  may  increase  Its  concentration  after  long  periods.  Further  work  to  analyze  and 
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heating  is  also  recommended.  In  addition,  experimental  data  required  for  using  the 
code  are  stated  and  some  experiments  are  suggested. 
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■  THIN  FILM  CONDUCTIVE  COATINGS  FOR  SURFACE  HEATING  AND  DECONTAMINATION 

1.  INTRODUCTION 

The  problem  of  contamination  of  solids  by  chemical  vapors  is  one  of  important 
consideration  for  the  design  of  military  vehic'e  components.  Typically,  a  vehicle  exposed 
to  undesirable  vapors  will  undergo  contamination  by  the  vapors  being  adsorbed  on  solid 
surfaces.  After  long  periods  of  exposure,  absorption  of  the  vapor  into  the  solid  will  take 
place.  This  can  cause  deterioratior.  of  materials  such  as  plexiglass  windows  and  may 
result  in  the  entire  unit  being  temporarily  nonfunctional. 

\n  the  present  analysis  we  examine  the  process  of  decontamination  by  heating  the 
plexiglass  substrates  with  imbedded  electrically  conducting  layers.  The  application  of 
such  heating  elements  for  the  purpose  of  deicing  is  well  known.  However,  little  is  known 
about  its  overall  effectiveness  for  the  removal  of  adsorbed  and  absorbed  contaminants. 
While  the  heat  input  supplies  energy  to  the  contaminant  molecules  and  sets  them  into  a 
free  state  (gaseous  state),  it  al?o  increases  the  solubility  of  the  contaminant  in  the  solid 
substrate.  It  is  therefore  necessary  to  carry  out  a  detailed  mathematical  analysis  to 
determine  the  effect  of  heating  on  adsorbed  and  absorbed  contaminants. 

In  the  current  development  we  first  present  an  overview  of  adsorption/desorption 
kinetics  and  then  develop  a  one-dimensional  model  for  the  removal  of  physically 
adsorbed  contaminants.  The  mathematical  analysis  provides  information  as  to  what  data 
are  needed  to  predict  the  performance  of  such  a  decontamination  system.  At  the  same 
time  some  typical  cases  have  been  run  to  simulate  such  prediction.  Also  included  are 
cases  for  which  the  increased  solubility  due  to  heating  may  cause  the  contaminant  level 

to  increase.  '  •'  '  ’  r':  *•  1  c*.  -Vf>  {  ^  \  \  v~  ;  i  y(-  :  —  \  , 

I  '  ;  •  ’ 

Further  to  the  conclusion  drawn  from  this  analysis,  recommendations  are  made  for 

future  studies.  In  particular,  calculations  for  optimal  heating  levels  to  minimi*.’  material 


7 


damage  by  heat  are  needed.  In  addition,  the  analysis  of  multiple-element 
decontamination  systems  and  further  d  pment  to  handle  chemisorption  and  multi- 
component  contaminants  are  very  much  needed  in  order  to  establish  design  guidelines. 


2.  KINETIC  THEORY  Or  ADSORPTION  AND  DESORPTION 

2.1  Fundamentals 

When  the  molecules  of  a  gas  strike  another  solid  or  liquid  surface,  they  tend  to 
cling  to  the  surface  tor  a  certain  length  of  time  befcre  going  back  to  the  gas  phase. 
Such  behavior  can  lead  to  a  high  concentration  of  these  molecules  on  the  surface.  This 
phenomenon  is  known  as  adsorption.  Molecules  may  be  adsorbed  either  physically  or 
chemically.  Physical  adsorption  takes  place  due  to  the  Van  der  Waals  forces.  Chemical 
adsorption  (or  chemisorption)  occurs  due  to  an  exchange  or  sharing  of  electrons  or  ions 
between  the  striking  molecules  and  the  molecules  of  the  surface. 

The  striking  molecules  when  adsorbed  enter  a  bound  slate  thereby  releasing  some 
energy.  This  release  of  energy  is  known  as  the  heat  of  adsorption,  if  an  adsorbed 
molecule  is  supplied  with  this  energy,  it  will  desorb  from  the  surface  and  go  into  the 
gaseous  state.  In  the  case  of  physical  adsorption,  the  heat  of  adsorption  is  about  10 
kilocalories  per  mole.  For  chemisorption  it  may  be  as  high  as  100  kilocalories  per  mole. 
The  adsorption  time  may  be  of  the  order  of  microseconds  for  physical  adsorption.  In  the 
case  of  chemisorption,  it  is  usual ly  several  seconds  butm3y,  in  some  instances,  be  as 
high  as  several  years.  For  both  regimes  of  adsorption  these  times  are  usually  several 
orders  of  magnitude  larger  than  the  thermal  vibration  time  of  the  molecules  of  the 
adsorbent.  Therefore,  the  number  of  molecules,  o,  adsorbed  at  any  giv«.i  instant  is  given 

o  •  nt  (2.1) 

where  n  is  the  number  of  molecules  striking  a  unit  area  per  unit  time  and  r  is  the 
adsorption  time. 

The  number  flux  at  the  surface,  for  an  ideal  gas,  is  given  by  the  following  well- 
known  relationship  from  kinetic  theory: 
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(2ttMRT)1/2 

where  N  is  the  Avogadro  number,  p  is  the  gas  pressure,  M  is  the  molecular  weight  of  the 
gas,  R  is  the  universal  gas  constant  and  T  is  the  temperature. 

The  adsorption  time  T  is  related  to  the  gas  temperature  and  the  heat  of  adsorption 


t  =  teQ'RT 


where  t  is  the  oscillation  period  of  the  molecules  and  Q  is  tne  molar  heat  of 
adsorption/desorption.  If  we  assume  that  all  the  molecules  striking  the  surface  are 
adsorbed,  then  by  combining  (2. 1-2.3)  we  see  that 


.  a  =  te  (2-4) 

(2ttMRT),/2 

Therefore,  for  given  vapor  adsorbate  and  adsorbent  surface,  at  constant  temperature,  we 
have  the  following  isotherm: 

o  =  kop  (2.5) 

where  k  is  a  constant. 

O  . 


The  number  of  molecules  per  unit  area  undergoing  desorption  by  thermal  excitation 
is  given  by 


nv  *  voe 


where  v  is  a  factor  of  proportionality  with  units  of  sec'1.  At  equilibrium,  therefore, 
n  *  n,.  and  v  *  1/t,  so  that  a  «  nteQ/RT  as  stated  earlier. 


if  jli  the  molecules  striking  a  sur/ace  are  not  adsorbed  and  a  fraction  a  is  reflected, 
then  n  is  reduced  by  a  fsctor  (1-a).  To. accommodate  for  this  we  define  an  apparent  time 
constant  x'  *  (1-a^)  and  use  it  instead  of  t. 
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The  adsorption  isobar  may  be  identified  from  (2.4)  by  holding  p  constant.  This  gives 


kIeQ/RT 


Since  the  T,/2  dependence  is  rather  weak,  we  may  write 


a  =  k,eL 


Similarly  at  constant  o  we  may  identify  the  isostere  as 


p  *  k2e-u/K' 

By  taking  the  logarithm  of  each  side  of  equation  (2.9)  we  may  write 


inp  ■  -  +  B 

RT  a 


(2  10) 


where  Ba  -  In  k2  is  a  constant.  This  equation  is  similar  to  the  expression  for  the  vapor 
pressure  as  a  function  of  temperature. 


*"[Po<T>l  -  -RT  *Bo 

By  defining  x  -  p/p  and  combining  (2.10)  and  (2.11)  we  have 


(2  11) 


in  x  -  4n(p/P0) 


(2  1.) 


Differentiation  with  respect  to  T  leads  to 


ain  x  Q-oo 
3T  RT2 

where  we  have  the  dependence  of  Q  and  Q  on  T  is  taken  to  be  negligible. 


(2.13) 


2.2  Langmiur's  Adsorption  Isotherm 

It  was  assumed  by  Langmuir  ^  that  adsorption  in  the  form  of  a  monomolecular 
layer  takes  place  at  a  constant  heat  of  adsorption  Q#.  It  was  further  assumed  that  the 
molecules  that  strike  the  covered  portiqn  of  the  interface  would  be  reflected 
instantaneously.  In  other  words,  the  residence  time  for  these  molecules  is  zero.  The 
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number  of  adsorbed  molecules  in  such  a  case  is  proportional  to  the  unused  sites  and  is 
given  by 

o  =  n't  =n(1-o/oQ)T  (2.14) 

where  oo  is  the  maximum  number  of  available  adsorption  sites  per  unit  area.  By  defining 
a  fraction  8  *  o/oo,  and  noting  that  from  (2.1)  and  (2.5)  we  have  nr  *  kop,  we  obtain 


1+k3p 


where  k,  =•  k  /cr  . 

J  O  O 

This  model  for  adsorption  is  suitable  for  low  pressures  and  is  particularly  applicable 
to  chemisorption  for  which  we  have  a  monomolecular  layer. 

2.3  Multilayer  Adsorption 

A  model  >*)r  multilayer  adsorption  was  suggested  by  Brunauer.  Emmett  &  Teller  2 
in  1938.  A  schematic  for  this  model  is  shown  in  Fig.  2.1  below.  Hero  the  regular  stacking 
of  the  molecules  is  shown  for  illustrative  purposes  only.  Actually,  the  molecules  would 
be  randomly  distributed  over  the  entire  surface  under  consideration. 


In  the,  model,  the  fraction  of  the  surface  covered  by  i  layers  of  molecules  is  denoted 
by  0(.  It  is  assumed  that  for  i  >  1.  the  heat  of  adsorption  Qa  is  constant  and  is 

identified  as  the  latent  heat  of  condensation  This  is  due  to  the  fact  that  the  energy 
exchange  ta>:es  place  between  like  molecules  as  it  does  in  the  case  of  condensation.  For 
i  -  0  the  heat  of  adsorption  while  capturing  a  monolayer  of  molecules  is  denoted  by  Qa. 
and  taken  to  be  a  constant  Generally  we  have  Q#  >  Q0 

In  order  to  maintain  the  fraction  as  constant,  the  number  of  molecules  adsorbed 
on  the  uncovered  portion  must  be  equal  to  the  number  of  molecules  desorbed  from  9, 
per  unit  time  per  unit  area.  This  would  correspond  to  an  equilibrium  state  for  Hence 
we  have 
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(2.21) 


Fig.  2.2 

A  schematic  of  total  heat  of  adsorption  calculation 


0,(0)  -  l(Q,  -  Q0)  >Qo01oo/N 
At  equilibrium  we  have 


9, 


1-x 

1-x  ♦  kx 


and 


1  -  9  -  -  0(1“*) 

°  1-x  ♦  lot 


With  a  little  algebra  it  can  be  shown  that 


0,(9)  •  {Q,  -  (0,  -Qo)*)0Oo/N  -  Q0Oo/N 
where  we  define  0  as 
0  -  0,  -  (0,  -  Q0)X 


(2.33) 


(2.34) 


(2.35) 


(2.36) 


(237) 


2.S  Limits  of  Validity  and  Further  Approximations 

For  x  <<  1  we  have  Q  ■»  Qa  as  expected  for  monolayer  adsorption.  As  x  1. 
Q  Qq  representing  a  near-saturation  condition. 


The  adsorption  isotherm  of  the  form  9  *  p  is  valid  only  for  extremely  low  partial 
pressures.  This  model  predicts  no  saturation  and  is  not  valid  for  p  -  p  .  The  Langmuir 
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k3P 


isotherm  of  the  form  ^  «  p  (or  9  *  (  ^  J  is  valid  only  for  monolayer  adsorption.  For 
Qa  >>  Q0  and  t,  <<  Tc  (i  >  1).  the  multilayer  results  reduce  to  the  Langmuir  model. 
Iri  the  case  of  pure  chemisorption  or  for  physical  adsorption  at  low  pressures,  the 
behavior  as  predicted  by  the  Langmuir  model  is  observed. 


2 

The  multilayer  adsorption  of  Brunauer  et  al.  is  based  on  constant  adsorption 
times  for  layers  corresponding  to  i  >  2  (above  the  monolayer).  Some  simplifications  are 
also  proposed  by  identifying  the  heat  of  adsorption  (which  is  taken  as  a  constant)  for 
higher  layers  witn  the  latent  heat  of  condensation  Q0.  The  identification  of  q  as  p0(T)  in 
(2  29)  is  an  additional  simplification  which  seems  quite  attractive.  However,  there  are 
situations  in  which  this  approximation  is  not  valid.  This  multilayer  adsorption  model  is 
applicable  over  a  wide  range  of  pressures. 

The  fundamental  principle  behind  the  decontamination  lies  in  behavior  of  9  with  . 
temperature.  Froai  (3  31)  it  is  clear  that  for  small  x,  9  is  nearly  proportional  to  x. 
However  po(T)  increases  with  T  and  hence  x  ■  p/p0  decreases  with  T.  Therefore  9 
decreases  with  T.  indicating  a  reduction  in  the  amount  of  contaminant  adsorbed  in  the 
surface  of  the  solid.  It  is  clear  that  by  raising  the  temperature  a  surface  may  be 
decontaminated. 

In  the  next  section  we  formulate  the  governing  differential  equations  for  substrates 
electrically  heated  by  imbedded  conductive  layers 
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3.  DECONTAMINATION  OF  SUBSTRATES  BY  ELECTRICAL  HEATING: 
FORMULATION 

3.1  Description  of  Problem 

Since  with  a,  rise  in  the  temperature  of  the  substrate,  decontamination  of  the 
surface  takes  place,  the  possibility  of  heating  by  imbedded  elements  (such  as  in  an 
automobile  windshield)  is  considered  here.  The  application  to  defogging  and  deicing  is 
well  known. 3»4t5,6  We  a(j0pt  here  a  one-dimensional  model  in  which  the  heat  and 
mass  flow  in  the  direction  parallel  to  the  plane  of  the  substrate  are  considered  negligible. 
At  this  point  we  can  define  a  specific  one-dimensional  time-dependent  problem. 

Let  us  consider  a  long  slab  of  thickness  L  which  is  exposed  to  an  environment 
containing  a  chemical  vapor  and  some  inert  gases.  An  electrically  heated  layer  of 
thickness  6  is  imbedded  in  the  slab  to  a  depth  h  (see  Fig.  3.1).  The  substrate  is  referred 
to  as  phase  1  and  the  heating  element  as  phase  2.  The  chemical  vapor  deposits  itself  on 
the  surface  of  the  substrate  by  adsorption  and  then  diffuses  into  the  bulk  of  the 
substrate. 

The  physico-chemical  processes  involved  are  as  follows: 

Mass  Transfer 

1.  Diffusion  and  convection  of  contaminant  vapor  in  the  environment  takes  place 
due  tc  the  wind  pattern,  or  to  the  motion  of  the  substrate.  This  sets  up  a 
velocity  profile  near  the  surface  of  the  substrate.  As  a  result,  convective 
transport  of  the  contaminant  to  or  from  the  surface  takes  place. 

2  Adsorption/desorption  of  the  vapor  at  the  solid  surfaces. 

3.  Diffusion  of  vapor  into  the  solid  phases  (1  and  2).  Chemical  reactions  within 
the  solid  are  not  being  considered. 


Solid  Phase  1 


Solid  Phase  2  (Conducting  Layer) 


Solid  Phase  1 


GAS  +  VAPOR  ("INSIDE") 

Fig.  3.1 

One-dimensional  model  of  substrate  with 
imbedded  heating  element 

Heat  Transfer 

1.  Heat  release  by  electrical  heating  within  the  conductive  layer. 

2.  Heat  conduction  in  the  solid  phases  1  and  2. 

3.  Heat  associated  with  adsorption/desorption  at  the  surface. 

4.  Thermal  diffusion  and  convection  in  the  gaseous  environment. 

Assumptions ' 

1.  Fluid  flow  processes  are  considered  only  for  their  effects  on  heat/mass 
transfer.  Simple  lumped  parameter  models  (i.e.,  heat  and  mass  transfer 
coefficients)  are  to  be  used  for  these  processes. 

2.  The  adsorbate  (solid)  is  infinite  in  length.  A  one-dimensional  formulation  is 
used  for  the  solid  phase. 

3.  A  uniform  volumetric  heat  generation  ratei  q ' "  is  considered  in  the 
conductive  layer  of  thickness  6.  Any  chemical  reactions  within  the  solid 
phase  and  the  latent  heat  release  thereof  is  not  considered. 

4.  The  concentration  of  the  sorbed  species  within  the  solid  is  linearly  related  to 
the  surface  concentration  of  the  adsorbed  layer  within  certain  limits. 
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5.  The  solid  surfaces  are  taken  to  have  reached  an  adsorption  equilibrium  with 
the  surrounding  gas-vapor  mixture.  At  this  equilibrium,  the  "fraction"  of  the 
surface  covered  (6)  depends  on  the  surface  temperature  Ts  and  the  partial 
pressure  pvs  of  the  vapor  adjacent  to  the  surface.  Furthermore,  at  equilibrium, 
pvs  is  purely  a  function  of  the  surface  temperature  and  the  heat  of  adsorption. 

If  pvs  is  different  from  psC0,  (the  far-field  partial  pressure)  then  vapor  transport 
takes  place  in  the  gaseous  phase  as  well. 

6.  For  tne  present  model,  only  physical  adsorption  is  treated  for  a  multi- 
molecular  adsorbed  layer.  The  heat  of  adsorption  for  the  first  layer,  Qa,  taken 
to  be  different  from  the  subsequent  layers,  for  which  it  is  taken  to  be  the 
latent  heat  of  condensation,  Q0.  For  the  first  layer  Qa  actually  varies  with  the 
fraction  covered  but  here  we  take  it  to  be  the  average  value. 

7  The  dependence  of  the  partial  pressure  with  temperature  is  given  by 


where  x  =  mvj  is  the  mole  fraction  of  the  vapor  in  the  air  adjacent  to  the 
surface. 

3.2  Governing  Equations 

At  this  point  the  problem  may  be  precisely  cast  into  a  mathematical  form  as  a 
closed  set  of  partial  differential  equations.  With  the  assumptions  made  in  §3.1,  the 
equations  for  heat  and  mass  balance  are  as  follows: 

Gas  Phase 

Heat  transfer  between  solid  surface  and  air: 

%  m  hg.(Ts  -  (3.2) 

Vapor  mass  flux  between  solid  surface  and  v/r. 

<g.v  -  hgmP9(m,.,  -mvV  (3  3) 


where 


with 


Nu  =  Nusselt  number 

kg  =  thermal  conductivity  of  gas 

niv  =  mole  fraction  of  the  vapor 

Dgi,  =  binary  diffusion  coefficient  between  vapor  and  air 
p  =  density. 


Equations  (3.2)  and  (3.3)  are  based  on  lumped  parameter  modeling  of  the  convective 
transport.  The  Nusselt  numbers,  Nut  and  Num  depend  on  the  external  flow  conditions. 
The  general  characteristics  of  these  Nusselt  numbers  are  available  in  most  heat  transfer 

textbooks  (see,  e.g.,  Burmeister  ^  ). 

Solid-Ges  Interface 


Heat  transfer  at  the  interface: 

%  " 
where 


Mass  transfer  at  the  interface: 
i  ♦  j  »  j 
where 

ia  *  mass  rate  of  adsorption  (g/cm2-sec) 
jsv  *  mass  flux  from  the  solid 


(3.7) 


(3  8) 


(3.9) 


3civ 

*+D - 

,v  3V 


y=-L 


and  jg  v  is  given  by  (3.3). 
Adsorption  Equilibrium 


From  (3.31) 


kx 

9  =  - 

(1-x)(1-x  +kx) 

where 


(3.10) 


(Q  -Q  )/RT 
k  =  e  3  0 


and  following  (2.37)  it  can  be  seen  that 


Q  -  Qo  =  (Qa  -  Q0)(1-x) 


(3.11) 


(3.12) 


By  employing  the  quasi-equilibrium  assumption  we  may  write  the  adsorption  heat  flux  as 

°o- 

qs*N9Q  (3.13) 

and 

°0-  . 

ia  *  N  0xMy  -moo9  (3.14) 

where  N  is  the  Avogadro  number,  a0  is  the  maximum  number  of  available  sites  per  unit 
area,  Mv  is  the  molecular  weight  of  vapor,  and  m  is  the  mass  of  one  molecule  of  the 
vapor. 

Solid  Phase  1 


Heat  transfer: 

1*T|_5ZT, 

a,  3t  3y? 

where  a,  is  the  thermal  diffusivity  and  T,  is  the  temperature. 


(3.15) 
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r 


(345) 


gives 


The  adsorption  equilibrium  given  by  (3.10)  and  (3.12).  when  non-dimensionalized. 


lex 


(1-x)(1-x  +kx) 


with 


(Q  -  Q  )/T 

k  »  e  *  »  * 


The  remaining  equations  may  be  written  as 


a  * 


3T1. 

3T, 

3t‘ 

aV*2 

■ 

K 

i 

at’ 

Le, 

3V*2 

32t2 

*  •  a 

at*  ’ 

♦<*• 

3v‘2  *  Q 

3C2. 

i  »s 

at’ 

♦o« 

le,  3V ** 

with  boundary  conditions 


c,w  -  c0  at  y  0.-1 


.  .*Ti 


3T, 


T1  -Tr3y.  •  J  y.  «yl»,-(h  ^ 


and 


*  *  ^C1»  ^C2y  ,  .  , 

and  initial  conditions 


T,  -  T2  -1 
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(3.46) 

(347) 

(348) 

(3.49) 

(3.50) 

(3.51) 

(3.52) 

(3.53) 

(3.54) 

\ 

\ 

\ 

(3.55) 


at  t  •  0 


We  next  examine  the  magnitudes  of  the  various  dimensionless  parameters  so  as  to 


identify  the  relatively  important  transport  mechanisms. 
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4.  IDENTIFICATION  OF  THE  IMPORTANT  TRANSPORT  MECHANISMS 

In  this  section  we  first  state  the  approximate  range  of  physical  properties  and 
external  conditions  for  typical  situations  of  practical  interest.  This  is  followed  by  an 
estimate  of  the  magnitudes  of  the  dimensionless  groups. 

4.1  Physical  Properties 

Plexiglass  substrate  (Poly  Methyl  Methacrylate) 


Property 

• 

Value 
at  378K 

Range 

Density, 

pi 

1.19 

1.0 

g/cc 

Specific  heat 

cpi 

0.37 

.5 

cal  'g-K 

Thermal  conductivity 

ki 

5x10'4 

1-IOxlO4 

cal/cm-sec-K 

Thermal  conductivity 

1.14x10'3 

1-IOxlO'3 

cm2/sec 

Mass  diffusivity  (02) 

D.v 

1x10‘® 

10'®-10'8 

cm2/sec 

Air  at  20°C 

Thermal  conductivity: 

*9 

06267X10"4 

.  cal/cm-sec-K 

Thermal  diffusivity: 

02216 

cm2/sec 

Mass  Oiffusivity  (02.N2) 

0.2 

cm2/sec 

Electrically  Conducting  Layer  (Indium  Oxide  *  Stannic  Oxide) 


Density: 

p2 

6.3 

g/cc 

Specific  heat: 

CP4 

02 

cal/g-K 

Thermal  conductivity: 

k2 

00136 

cal/cm-sec-K 

Thermal  diffusivity: 


cm2sec 


a2  1.08  x  10~2 

Mass  diffusivity:  D2  Not  Available  (consider:  D2v/Dl¥  *  10‘\  1  and  10) 

Thickness  of  layer  5  typically  1000  A*  10Jcm,  (consider:  1  -  10  x  10'Jcm) 


Heats  of  Adsorption 

The  ambient  temperature  is  taken  to  be  T^  »  20°C  (293  K).  Assuming  Qa  »  3 
kcal/g-mole,  we  obtain: 

Q.  "Q./RT«  a  5 2  (4. i) 

With  this  as  an  approximate  estimate,  a  range  Q#  *  1-25  is  considered. 

Similarly,  with  Qo  *  1  kcal/g-mole,  Qo  -  1-5  is  considered.  . 

Heating  Levels 

The  maximum  heating  levels  quoted  in  the  literature  are  of  the  order  of  4 
cal/cm2-sec.  Therefore,  the  volumetric  heating  rating  (for  5  »  10'3cm)  is 


.«  qL2 

q  »  - 

We  consider  q  * 


»  60 

-  1-100 


(4.2) 


Solubility  Parameter 

The  solubility  parameter  e  is  estimated  as  follows: 

At  the  surface  we  have  c1v  •  €0.  As  9  takes  on  large  values  (say,  10)  the  solid 
phase  will  also  approach  saturation  (c^  ♦  1).  Therefore,  we  take  r  ■  1/10  •  0.1. 


Adsorption  Depth 


The  dimensionless  adsorption  depth,  H  -  moo/Lc0  is  estimated  by  examining  the 
solubility  of  N2  in  Poly  Ethyl  Methaclylate.  At  25°C  the  solubility  is 

s  1  7.5  x  10'2c.c.gasSTP/c.c.substrate  (4.3) 

where  s  is  defined  by  c1v  3  sp.  Assuming  a  partial  pressure  of  p  »  0.1  atm.  we  find  the 
volumetric  concentration  c,  to  be 

IV 

c1v  3  7.5  x  10~3c.c  (STP)/c.c.  (4.4) 

In  units  of  mass  concentration  this  is 

c1v  3  9  x  10'6g/c.c.  *  10‘5g/c.c.  (4.5) 

We  therefore  take  cq  a  10'5  g/c.c..  If  the  gas  is  more  soluble,  may  be  as  large  as  1<T3 
g/c.c..  For  very  low  solubilities  it  may  be  10*7  or  10~8  g/c.c. 

4.2  Dimensionless  Groups 

The  maximum  number  of  available  sites  is  of  the  order  of  a  ?•  10,4/cm2.  The 
parameter  H  is  therefore  given  by 

mo„ 

H  -  -  10'4  (4.6) 

0 

where  l  is  taken  to  be  1  cm,  and  m  3  4.67  x  10'23g  for  Nz. 


negligibly  small.  Furthermore,  due  to  the  very  small  mass  diffusivity  of  the  solid  (10~6  - 
10-7  cm2/sec)  compared  to  that  in  the  gaseous  phase  (0.1  cm2/sec),  the  mass  transfer 
relationship  betwen  the  surface  and  the  ambient  reduces  to  the  quasisteady  relation 

c1v  m  (4.8) 

where  9  is  the  surface  coverage  a t  the  adsorption  quasi-equilibium. 

The  above  linear  relationship  can  be  generalized  by  considering 
phenomena  within  the  solid  and  a  Langmuir  type  absorption  isotherm, , 

e0 
c  * 

1  ♦  60 

can  be  employed. 

The  results  of  the  estimates  of  these  dimensionless  groups  are  summarized  in  Table 
4.1. 


saturation 


(4.9) 
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Table  4-1.  Summary  of  estimate  of  magnitudes  of  dimensionless  parameters 


Dimensionless  Parameter 

Range 

Comments 

Nu, 

0(100) 

Included  in  the  model 

Num 

0(100) 

This  is  large,  but  irrelevant 
since  gas-phase  mass  transfer 
does  not  affect  problem 

A 

10'9 

Negligible 

Q 

O 

1-5 

Included 

Q 

a 

1-25 

Included 

•  •- 

q 

1-100 

Included 

H 

10'6  -  10’1 

Irrelevant 

e 

<1 

Strongly  temperature  dependent. ' 
Consider  values  10"4  -  1 

0.001 

Consider  0.001,  0.01,  0.1 

h 

0  <  h  <  1 

0.05  -  0.5 

-  ' 

♦o9 

105  -  10® 

Le, 

102  -  105 

♦as 

1  -  10 

K 

10  -  100 

♦os 

— 

No  data  available,  try  0.1.  1.  10 
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S.  RESULTS  AND  DISCUSSION 

The  non-dimensional  governing  differential  equations  together  with  the 
simplifications  discussed  in  §4  have  been  programmed  for  a  finite-difference  solution.  A 
listing  of  the  program  is  given  at  the  end  of  this  section.  The  input  parameters  required 
to  run  the  program  are  also  given  in  the  listing.  For  the  sake  of  clarity  these  parameters 

are  listed  below  in  the  mathematical  notation  used  in  §4. 

Program  Mathematical 

Notation  Notation 

DELTA  6“  *  5/L,  thickness  of  heating  layer 

H  h  3  h/L,  heater  location 

HEAT  q*  *  q'^/kT^  volumetric  heating  rate 

:V 

QA  Qa  =•  CWRT,,,,  heats  of  adsorption 

QO  Q  3  Q /RT_  heats  of  condensation 

O  o  w 

TINIT  Too/T«o  *  1,  initial  temperature  =?  ambient 

CINIT  0,  initial  concentration  in  the  bulk 

.  * 

PVBOT  p  .  *  p  , ,  dimensionless  partial  pressure  of 

y  CO  *  v  VL  k 

vapor  on  the  "inside"  (i.e..  bottom  side  of  substrate.) 

*  * 

PVTOP  p  *  p  ~  dimensionless  partial  pressure  of  vapor 

y,09  I  O  yu 

on  the  "outside*  (i.e..  top  side  of  substrate.) 

ANUSTT  Nut0,  thermal  Nusselt  number  on  the  "outside* 

ANUSTB  Nu(L,  thermal  Nusselt  number  on  the  "inside" 

FIKG  <>k9  *  kg/kr  thermal  conductivity  ratio  (gas/piexiglass) 

FIKS  $kt  *  k2/k,.  thermal  conductivity  ratio  (conductive  Isyer/plexiglass). 

FIALFS  $aj  *  a2/ar  thermal  diffusivity  ratio  conductive  layer/plexiglass) 

FIDG  $0g  »  Dg/D,,  mass  diffusivity  ratio  (gas/plexiglass) 

FIOS  $0f  *  D2/Dt,  mass  diffusivity  ratio  (conductive  layer/piexiglass) 
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TM'/VT'V  If  \m 


>.m  t  r»  r  »  yr 


\|i|*ii  ^  i  ^ 1 


T*T*rrrrw 


■■  !■•  ■  ?"«?"•  i  •.  v 


ALEW1  Le,  ■  o^/D^  Lewis  number 

EPS  £  =  moo/^c0,  solubility  parameter 

The  results  from  the  various  runs  of  the  program  have  been  plotted  in  Figs.  5. 1-5.5. 
Here  we  discuss  each  case  in  detail. 


In  Fig.  5  1.  the  temperature  profile  at  various  times  is  shown.  The  heater  is  placed 
at  y  ■  -0.25.  The  thermal  parameters  for  this  plot  are  q  *  50.0,  6  *  0.001,  $ks  *  50.0, 
Nut0  3  200.0  and  NutL  3  10.0.  The  plot  shows  the  following  important  features: 

1.  The  region  near  y  3  0  reaches  a  steady  state  faster  than  the  rest  of  the 
substrate.  This  is  because  of  the  higher  Nusselt  number  and  the  shorter 
distance  from  the  heating  element. 

2.  More  heat  leaves  through  the  surface  y  3  0  than  y  3  -1.  This  is  owing  to  the 
relatively  lower  thermal  resistance  of  the  region  -h*  <  y  <  0  than 

-1  <  y  <  -h  . 

3.  A  large  Nusselt  number  causes  the  corresponding  surface  to  be  cooler  and,  as 
a  result,  leads  to  higher  adsorption.  The  heating  is  therefore  wasteful.  If  the 
Nut  is  controllable,  then  it  should  be  minimized  so  that  very  high  heating 
levels  are  not  needed. 

4  The  maximum  steady  temperature  in  this  case  is  T  /T*  *  1.5.  Assuming 
Tg,  3  300K,  we  have  Tma>;  «■  450K.  At  such  high  temperatures  the  plexiglass 
would  deteriorate. 


In  Fig.  5.2  the  dimensionless  surface  coverage  of  the  contaminant  on.  the  outside, 
(0o)  and  on  the  inside  (0L)  are  shown  as  functions  of  time.  The  parameters  for  the  graph 
are: 


H 

S’1 

s  . 

•  A 

q 

Nu«„ 

Nua 

Run  1 

0.5 

0.0001 

10 

0.05 

50 

200 

10 

1 

Run  II 

0.25 

0.001 

50 

0.25 

50 

200 

10 

5 

Run  III 

0.5 

,  0.001 

50 

0.25 

50 

10 

10 

5 

Run  IV 

0.5 

0.001 

50 

0.25 

10 

2 

2 

5 

in  addition  we  use  Qa  3  5.0, 

CL  •  1.0, 
0 

m  0  5 

and  pwl  - 

0.01. 
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variation  of  T*  takes  place  at  the  heater  location 


Temperature  distribution  within  the  substrate  at  various  times 


Fig.  5.2 

Dimensionless  surface  coverage  at  various  times 


Here  we  observe  the  following  features: 


1.  Curves  I  and  II,  between  which  there  is  no  systematic  change  in  the  thermal 
parameters,  both  reach  the  same  steady  state  value  for  0Q  and  0L.  This  is 
because  the  groups 


1  +  Nu'°  V: 

1  i  +  -Nu(L  4>kg0-h’) 

and 


(5.1) 


<3 5  5  *ks 

*kg(Nuto  +  NuttQ1> 


have  approximately  the  same  values  for  these  curves. 


(5-2) 


2.  The  curve  II  corresponds  to  a  larger  thermal  diffusivity  than  curve  I.  It 
therefore  approaches  a  steady  state  faster.  For  curve  IV  the  Nusselt  number 
is  lower  that  curve  III.  The  eigenvalues  determining,  the  rate  of  thermal 
transport  are  smaller  for  curve  IV  and  the  transport  process  lasts  longer  in 
this  case. 

3.  The  higher  level  of  adsorption  on  the  outside  is  due  to  larger  partial  pressure 
of  the  vapor  on  the  outside. 


In  Fig.  5.3,  the  effects  of  partial  pressure  of  the  vapor  anti  the  heats  of  adsorption 
on  the  fraction  covered  are  shown.  The  plots  correspond  to  fixed  values  of  the 
parameters  Q,  and  Q2  (Q1  *  0.0317,  Q2  »  11.5556)  or  fixed  values  of  the  parameters 

T,  *  1  ♦  Q2  (5.3) 

and 

Tt  •  1  ♦  Q,Q2.  (5.4) 

1.  We  find  that  0L  <  90  because  Tu  >  Ta.  This  is  due  to  the  convective  cooling 
which  occurs  on  the  outside 

2  Each  curve  shows  a  steep  linear  portion  for  small  pressures,  a  flat  portion  for 
moderate  pressures,  and  saturation  (0  ♦  •)  as  the  partial  pressure  approaches 
saturation. 
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Diroens 


3.  At  large  Q  the  linear  and  the  flat  portions  are  separated  around  9-1.  this 
implies  that  for  large  Qa  (such  as  in  chemisorption)  a  monolayer  is  formed 

A 

first  until  9-1.  Subsequently,  more  layers  build  up.  For  Qa  Qo,  all  layers 
•nay  be  formed  simultaneously. 


4.  At  large  Qq.  the  saturation  phenomenon  is  delayed.  This  happens  because  the 

|Q  -Ql/T 

parameter  k  =  e  a  0  s,  which  signifies  the  ratio  of  adsorption  times 


between  the  first  and  the  subsequent  layers,  decreases. 


In  Fig.  5  4,  the  concentration  profile  within  the  solid  is  plotted.  The  solubility  has 
been  assumed  to  be  the  same  in  both  the  substrate  and  the  heating  element.  Also  the 
ratio  of  the  two  diffusivities  is  taken  to  be  unity. 


The  parameters  are:  L  *  0  5.  5  -  0.001,  *  5.0,  $Ds  *  1.0,  $ko  *  50.0,  le,  -  500, 

$kg  -  0.25,  $Dg  -  5  0  x  10s.  Qa  -  5.0,  Qo  *  1.0.  q  -  50  0.  Nu!o  -  10  0.  Nu,L  -  200  0.  e  -  0.1. 
Pvo  -  0.5  and  qvL  -  0  01 


The  plot  shows: 

1.  The  time  taken  to  reach  mass  transfer  steady  state  is  approximately  equal  to 
Le,  x  (time  for  thermal  steady  state). 

2  For  short  times,  diffusion  occurs  from  both  ands  and  the  effects  from  each 
end  grow  independently  until  they  interact. 

3.  As  time  increases,  the  concentration  profits  becomes  monotonic  and  a 
straight  line  is  obtained  for  constant  mass  diffusivity.  At  steady  state,  a 
steady  stream  of  vapor  diffuses  from  the  higher  concentration  side  to  the 
lower  concentration  side 

4.  The  surface  concentrations  are  given  by 


c0  “  £  V<1  *  at  V  "0 
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cl  *  £®l/^  *  at  v '  *1 


(5.6) 


Here  9o  and  0L  are  functions  of  the  thermal  parameters  Q,  and  Q2.  the  heats 
of  adsorption  Qa  and  Qq,  and  the  partial  pressures  pwo  and  p^L  Thus  after  the 
thermal  steady  state  is  reached,  cq  and  cL  remain  constant  due  to  thermal 
equilibrium. 

In  Fig.  5.5.  we  have  plotted  the  variation  in  the  steady  state  bulk  concentration 
(cbulk  *  J°jC  dy  )  as  a  function  of  the  heating  rate  q  Also  plotted  in  the  same  graph  is 
Tmjj  j  Tma„/TaB  as  a  function  of  q  .  The  important  feature  incorporated  here  is  that  the 
solubility  parameter  e  is  taken  to  be  temperature  dependent.  It  is  given  by 


Q  (1  - 1/T  i 

e  -  e0e  *  * 


(57) 


where  *  Q/RToo  is  the  heat  of  solution.  Since  the  solubility  changes  with  changing 
temperature,  the  heat  of  solution  plays  a  role  The  program  was  modified  to  include  this 
addition  parameter  The  plot  Corresponds  to  the  following  values  of  the  parameter: 


l  »  0  8. 

6  -  0001. 

-  50  0. 

♦*9  -  0  25. 

Nu  -  200  0. 

to 

Nu,l  -  10  0. 

-  0  1. 

O 

Q  -  10  0. 

• 

Q  -  10. 

O 

p  -0  5. 

VO 

p  .  -  001. 

vl 

o 

i 

o 

o 

¥i 

a 

The  results  exhibit  the  following  features 

1.  For  very  low  heat  of  solution.  Q.  the  bulk  concentration  decreases 
monotonically  with  the  heating  level  q  It  may  be  noted  that  the  duration  of 
heating  does  not  determine  the  level  of  contamination  after  the  attainment  of 
thermal  and  mass  transfer  steady  states  For  this  reason  the  level  of  initial 
contamination  or  the  initial  temperature  do  not  affect  cbu(fc. 

2  For  a  reasonably  large  heat  of  solution,  while  the  surface  coverage  (0O  and  0L) 
decreases  with  increasing  temperature,  the  dissolved  contaminant  in  the  bulk 
increases.  This  usually  happens  when  the  heating  levels  are  low  and  the 
change  of  0Q  and  0L  are  not  as  rapid  as  that  of  the  solubility.  If  the  outside 
environment  is  very  cold  and  if  the  convective  cooling  is  strong,  the  interior 
of  the  substrate  may  be  very  hot.  but  the  surface  will  remain  fairly  cool.  As  a 
result,  the  heating  may  hot  substantially  remove  the  surface  contaminant  and 


'r 1  .»  r  J”  „■ 1  r 1  l"  tj  l*  «■'  ■  v.  .r  rr  v 


at  the  same  time  it  will  increase  the  solubility.  This  will  lead  to  increased 
contamination  if  heating  is  sustained  for  long  periods. 

3.  An  important  consideration  for  design  would  be  the  maximum  temperature 
reached: 

Tc  +TLh  /(1-h  )  +q  6  h  <J»ks 
Tmax  1  +h  /(1-h  ) 

If  =  300K,  then  for  a  material  such  as  plexiglass,  we  would  require 
T  <  1  -5  t0  avoiding  thermal  damage. 


(57) 
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6.  RECOMMENDATIONS  AND  CONCLUDING  REMARKS 

Here  we  further  discuss  the  feasibility  of  decontamination  by  heating  and  the 
limitation  of  the  effectiveness.  We  also  make  recommendations  for  experimental 
measurements  for  obtaining  the  required  data  in  order  to  run  the  computer  code.  In 
addition,  we  make  several  recommendations  so  that  all  of  the  important  physical 
phenomena  may  be  incorporated  into  the  computer  program,  and  that  concrete  design 
criteria  may  be  extablished  for  engineering  use. 


6.1  Recommendations  for  Experiments 

The  required  dimensionless  parameters  are  all  given  in  S5.  Most  of  the  data 
heeded  to  calculate  these  parameters  are  available  except  for  some  information  on 
adsorption/desorption  kinetics  and  thermal  permeability  of  the  solid  substrate.  In 
particular,  the  following  measurements  are  required: 

a.  Adsorption  Isotherm:  The  functional  dependence  of  the  heat  of  adsorption  on 
the  temperature  [i.e.,  Q(T)]  is  a  necessary  parameter.  For  this  measurement 
we  may  take  a  solid  sample  and  heat  it  up  in  order  to  remove  all 
contaminants.  Then  at  a  given  partial  pressure  of  the  contaminant,  expose  it 
for  a  short  period  of  time  (say,  2-3  minutes)  while  holding  its  temperature 
fixed  at  T.  By  subsequently  reheating  the  sample  and  measuring  the  volume  of 
gaseous  contaminant  that  desorbs,  the  surface  coverage  6  may  be  calculated. 

Then  we  could  calculate  the  constant  using  x  -  p/pQ.  By  using  (2.32)  and 
noting  that  Qa  is  approximated  to  be  the  latent  hest,  of  condensation,  Qa  can 
be  found.  Equation  (2.37)  then  yields  the  value  of  Q. 

b.  Permeability:  The  mass  diffusion  coefficient,  D1v  and  the  maximum  solubility 
c0  are  the  other  physical  properties  that  have  to  be  measured  at  various 
temperatures.  For  this  measurement  one  should  take  a  thin  sheet  of  the 
uncontaminated  (or  decontaminated)  substrate  and  expose  it  to  a  contaminant 
on  one  side.  By  letting  the  contaminant  penetrate  to  the  other  side  and 
measuring  the  concentration  profile  as  a  function  of  time  one  may  calculate 
the  difusion  coefficient.  When  the  substrate  saturates,  the  maximum 
concentration  cQ  may  be  calculated. 
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6.2  Recommendations  for  further  Work 


The  present  analysis  is  based  on  one-dimensional  heat  and  mass  diffusion  in  the 
}  substrate  having  constant  physical  properties.  Only  a  simple  variation  of  the  solubility 

with  temperature  has  been  treated  (Fig.  5.5)  for  the  purpose  of  examination  of  adverse 
effects  of  heating.  Since  many  of  the  other  properties  are  held  constant  the  accuracy  is 

i 

I  *  likely  to  be  poor  In  order  to  establish  design  criteria  for  quantities  such  as  optimal 

heating  rate  or  for  optimal  heater  location  further  work  is  clearly  necessary  We 
recommend  the  following: 

i.  Complete  Variable  Property  Analysis:  The  mass  diffusivity.  vapor  solubility  and 
the  heats  of  adsorption  have  to  be  considered  as  functions  of  temperture.  As 
noted  earlier,  property  variations  can  cause  the  decontamination  system  to  be 
not  only  ineffective,  but  counterproductive. 

i 

ii.  Multi-Element  Optimal  Heating:  For  thick  substrates,  a  great  deal  of  heating 
may  be  necessary  in  Order  to  achieve  the  required  surface  temperature.  This 
can  cause  heat  damage  as  well  as  contamination  due  to  increased  solubility. 

The  solution  to  the  problem  includes  the  possibility  of  having  two  or  more 
elements  in  order  to  have  a  more  even  distribution  of  heat.  Optimal 
|  positioning  in  this  case  would  be  quite  important  especially  if  the  desorption 

characteristics  are  different  on  each  of  the  two  sides  of  thn  substrate. 

iii.  Multicomponent  Adsorption:  The  current  theory  only  deals  with  a  single 
component  of  contaminant  vapor.  In  real  situations  two  or  more  different 
species  of  vapor  may  be  present.  A  generalization  to  account  for  this  is 
needed. 

iv.  Chemisorption  and  Chemical  Reactions:  The  model  developed  so  far  only 
takes  care  of  physical  adsorption.  Various  modifications  are  necessary  in 
order  to  deal  with  chemical  adsorption.  Also,  there  may  be  some  important 
effects  due  to  chemical  reaction  in  the  bulk  of  the  substrate 

v.  Two-Oimensional  Analysis:  Our  one-dimensional  analysis  utilizes  constant 
average  transport  parameters  in  the  gaseous  phase  (Nu().  For  practical 
circumstances  the  momentum  and  the  thermal  boundary  layer  thicknesses 
vary  along  the  surface  of  the  substrate.  Further  modification  to  allow  for  this 
is  required 


6.3  Conclusions 


With  a  complete  analysis  as  recommended  accurate  design  criteria  can  be 
established  on  a  sound  scientific  basis.  A  high  degree  of  accuracy  may  be  necessary  for 
some  critical  cases  for  which  heating  levels  may  have  to  be  pushed  to  the  maximum 
limits  As  noted  in  items  (i)  and  (ii)  careful  design  is  needed  for  substrates  exposed  to 
certain  extreme  conditions  under  which  heating  may  actually  increase  the  contaminant 
level,  it  is  also  clear  from  this  study  that  heating  can  be  used  to  preclude  a  substrate 
from  being  contaminated.  Here  again  the  optimal  heating  levels  for  various  external 
circumstances  need  to  be  established  on  the  basis  of  practical  considerations. 

in  closing,  we  affirm  that  there  is  definite  merit  in  conductive  film  heating  for 
decontamination  or  for  precluding  contamination.  With  the  acquisition  of  recommended 
data  and  with  further  generalization  of  the  computer  code,  workable  systems  may  be 
designed  with  the  least  possible  experimentation. 
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GLOSSARY 


concentration  of  absorbed  vapor  within  the  solid 

maximum  vapor  concentration  that  is  dissolvable  in  the 
solid  (This  is  related  to  the  available  absorption  sites 
per  unit  volume) 

specific  heat 

mass  diffusion  of  vapor  (in  ambient  gas  or  solid) 
heater  depth  from  outside  surface 
dimensionless  adsorption  depth  [equation  (3.28)] 
thermal  conductivity 
thickness  of  the  slab 
Lewis  number  (  =  o/D) 

mole  friction  of  the  contaminant  vapor  molecules 

Nusselt  number 

partial  pressure  of  vapor 

volumetric  heating  rate 

molar  heat  of  adsorption  for  the  first  layer  of  adsorbed  molecules 

molar  heat  of  condensation 

molar  heat  of  solution 

universal  gas  constant 

time 


temperature 


• ,  .  ■  \ 

/  '  '  - 

Ts  -  surface  temperature 

To,  -  ambient  temperature 

x  -  dimensionless  partial  pressure  (also  equal  to  mote  fraction) 

y  -  coordinate  normal  to  the  slab 

Greek  Letters 

a  -  thermal  diffusivity 

6  -  thickness  of  heating  element 

A  -  coupling  parameter  governing  the  thermal  transients  of  the 

slab  surface  [Equation  (3.27)} 

e  -  solubility  parameter  which  links,  surface  coverage  and  bulk 

concentration  at  the  surface  [c  *  £0/(1  £0)1 

♦  -  property  ratio 

5  -  density 

a  -  dimensional  surface  coverage 

aa  -  available  i  adsorption  sites  per  unit  area 

9  -  dimensionless  surface  coverage  (  ■  Q/oJ 

Subscripts 

0  -  upper  sui-face 

1  *  substrate  material  (plexiglass) 

2  -  heater  material  (Indium- Tin  Oxide) 

•  -  Ambient  Far-stream  quantity 

0  -  mass  diffpsivity  ratio 

g  -  for  prepefrty  of  gas-vapor  mixture 

k  -  thermal  condictivity  ratio 
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lower  surface 


mass  transfer 

at  the  slab  surface;  property  ratio  between  solids  2  and  1 
thermal 

contaminant  vapor 


tnermal  diffusivity  ratio 
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APPEN01X 


THIS  PROGRAM  SOLVES  FOR  THE  DIFFUSION  OF  HEAT  AND  CONTAMINANT 
VAPOR  WITHIN  A  SOLID-SUBSTRATE  (PLEXI-GLASS)  WITH  AN  EMBEDDED 
CONDUCTIVE  LAYER  (  INDIUM-TIN  OXIDE).  THE  NUMERICAL  PROCEDURE 
IS  BASED  ON  THE  IMPLICIT  FINITE  DIFFERENCE  PROCEDURE.  THE 
PROGRAM  IS  CURRENTLY  SET  FOR  A  MAXIMUM  OF  1000  INTERVALS  FOR 
SPACIAL  DIFFERENCING.  IF  'N'  IS  THE  NUMBER  OF  INTERVALS  REQUIRED 
THE  ARRAYS  HAVE  TO  BE  DIMENSIONED  AS  'N+l'. 


THE  INPUTS  TO  THE  PROGRAM  AREJ 


NR»NW  -  DEVICE  NUMBERS  FOR  READING  DATA  AND  WRITING 
OUTPUT.  THESE  ARE  SPECIFIED  THROUGH  THE  DATA 
STATEMENT  IN  THE  BEGINNING  OF  THE  PROGRAM. 
MAXSTP  -  MAXIMUM  NUMBER  OF  TIME  STEPS 
IL 1ST  -  NUMBER  OF  TIME  STEPS  AFTER  WHICH  SURFACE 

TEMPERATURE  AND  CONCENTRATION  AND  INTEGRATED 
BULK  CONCENTRATION  WILL  BE  PRINTED. 


IGRID1  -  NUMBER  OF  GRID  SPACES  BETWEEN  THE  HEATER  AND 
INSIDE  SURFACE 

IGRID2  -  NUMBER  OF  GRID  SPACES  IN  THE  HEATING  LAYER 
IGRID3  -  NUMBER  OF  GRID  SPACES  BETWEEN  THE  OUTSIDE 
SURFACE  AND  THE  HEATING  LAYER 
IPROF  -  NUMBER  OF  TIME  STEPS  AFTER  WHICH  TEMPERATURE 
AND  CONCENTRATION  PROFILES  ARE  REQUIRED 
DELTA  -  THICKNESS  OF  HEATING  LAYER 
H  -  DEPTH  AT  WHICH  HEATER  IS  LOCATED 

DT  -  TIME  STEP »  FOR  STUDY  OF  THERMAL  TRANSIENTS  A 

TIME  STEP  OF  1/MAXSTP  IS  RECOMMENDED.  FOR  THE 
STUDY  OF  MASS  TRANSIENTS  WITHIN  THE  SOLID*  A 
TIME  STEP  OF  (LEWIS  #)/MAXSTP  IS  RECOMMENDED. 
HEAT  -  DIMENSIONLESS* VOLUMETRIC  HEATING  LEVEL 
OA  -  DIMENSIONLESS  HEAT  OF  ADSORPTION 

QO  -  DIMENSIONLESS  HEAT  OF  CONDENSATION 

TINIT  -  DIMENSIONLESS  INITIAL  TEMPERATURE 

CINIT  -  DIMENSIONLESS  INITIAL  CONCENTRATION 

PVBOT  -  DIMENSIONLESS  PARTIAL  PRESSURE  OF  VAPOR  ON 

THE  INSIDE 

PVTOP  -  DIMENSIONLESS  PARTIAL  PRESSURE  OF  VAPOR  ON 
THE  OUTSIDE 

ANUSTT  -  THERMAL  NUSSELT  NUMBER  ON  OUTSIDE 
ANUSTB  -  THERMAL  NUSSELT  NUMBER  ON  INSIDE 


FIKG  -  THERMAL  CONDUCTIVITY  RATIO  (GAS/PL.  GLASS) 

F IKS  -  THERMAL  CONDUCTIVITY  RATIO  (COND.  LAYER/PL.  GLASS) 

FIALFS  -  THERMAL  DIFFUSIVITY  RATIO  (COND.  LAYER/PL . GLASS ) 
FIDO  -  MASS  DIFFUSIVITY  RATI0(GAS/PL.  GLASS) 

FIDS  -  MASS  DIFFUSIVITY  RATIO(COND.  LAYER/PL.  GLASS) 

ALEW1  -  IEWIS  NUMBER  OF  PLEXI-GLASS 

EPS  -  SOLUBILITY  PARAMETER 
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c* 

c*  THE  GIVEN  INPUT  DATA  IS  PRINTED  OUT.  HERE  'NW'  REFERS 

C*  TO  THE  OUTPUT  DEVICE  NUMBER. 

C* 

WRITE ( NU  » 300  >  IMAX »H » DELTA. DT 

WRITE (NW. 301)  FIALFS.FIDS.FIKS.ALEW1 .FIKG.FIDG 
WRI TE ( NW » 302 )  OA. 00. HEAT 
WRITE (NW.303)  ANUSTB » ANUSTT 
WRITECNW.304)  EPS 

WRITE ( NW . 305 )  TINIT.CINIT .PVBOT.PVTOP 

300  FORMAT  (10X.  'NUMBER  OF  GRID  PTS  =  ' » 15/ .10X.  'DEPTH  OF  HEATER  =  ', 
.  *F8.5/» 10X. 'THICKNESS  OF  HEATED  LAYER  *  '.F8.5/.10X. 

♦'TIME  STEP  =  ' »F8«5//) 

301  FORMATdOX. 'THERMAL  DIFFUSIVITY  RATIO  (S2/S1)  =  '  »F10.5/»  10X. 

* ' MASS  DIFFUSIVITY  RATIO  (S2/S1)  =  '.F10.5/.10X. 'THERMAL  ' 

♦.'CONDUCTIVITY  RATIO  (S2/S1)  =  ' .F10.5/.10X. 'LEWIS  NUMBER 

*  OF  SOLID  1  =  ' ,D12. 4/. 10X. 'THERMAL  CONDUCTIVITY  RATIO  <G/S1>=  ' 
♦.D12.4/.10X. 'MASS  DIFFUSIVITY  RATIO  (G/Sl)  =  '.D12.4//) 

302  FORMAT (10X. 'HEAT  OF  ADSORPTION  =  '  »F10.5» 10X. 'HEAT  OF' 
♦.'CONDENSATION  =  ', F 10. 5/ . 10X »' HEATING  RATE  =  '.F10.5//) 

303  FORMAT ( 10X. 'NUSSELT  NUMBER  (THERMAL)  ON  INSIDE  =  '.F10.5/ 

♦ . 10X . 'NUSSELT  NUMBER  (THERMAL)  ON  OUTSIDE  =  '.F10.5//> 

304  FORMAKl OX. 'PARAMETER  EPSILON  =  '.D12.4//) 

305  rORMATdOX, 'INITIAL  TEMP  =  ' .F10. 5. 10X. ' INITIAL  CONC  =  '.F12.8/ 
♦10X. 'VAPOR  FRACTION  INSIDE  = ' .F10.5. 10X, 'VAPOR  FRACTION'. 

♦'  OUTSIDE  »  ' .F10.5// ) 

C« 

C**  INITIAL  CONDITIONS  ARE  SPECIFIED. 

C** 

DO  1  1*1. *  IMAX 
OTH(I)  =  0.0 
0D( I )  a  o.O 
Cd)  =  CINIT 
t  KI)  *  TINIT 

DO  2  I=IHLOW.IHUP 
QTHd)  *  HEAT 

C«4 

C***  THE  COEFFICIENTS  FOR  THE  TRI-DIAGONAL  FINITE  DIFFERENCE 

C***  EQUATIONS  ARE  CALCULATED. 

C*** 

F1ALF1  *  1.0 
FIALF2  *  FIALFS 
FID1  *  1 . 0/ALEW1 
FID2  =  F IDS/ALEW1 
Cl  *  DT/(DYltDYl ) 

C2  *  DT/(DY2tDY2> 

C3  ~  DT/(DY34DY3> 

DM1  *  C 1 ♦F I D1 
DM2  *  C24FID2 
MM3  ==  C34F I D1 
F  IMK  *  FIDS 

GMH  *  AMB/FID1  ♦  AMB$AMB4F I MK/F I D2 

GMT  -  AMT/FID1  F  AMT*AMT ♦FIMK/FID2 

BT1  *  C 1 ♦F IALF 1 

DT2  *  C24FIALF2 

DT3  *  C3tFIALFl 

FI  TK  FIKS  , 

GIB  *  AMD/FI ALT i  F  AMB^AMBfF I TK/F I ALF2 
GTT  *  AM 7 /Fi ALU  +  AMT$AMT*F I Tk/F I ALF2 
CALL  C0£FM(BM1  »BM2.BM3»FIMK»GMB.GMT.SUBM»DIAH.SUpM) 

CALL.  , CQEFM( BT 1 » DT2» BT3.F ITK.GTB.GTT.SUBT.DIAT.SUPT > 


IMPLICIT  REAL*8(A-HfO-Z> 

DIMENSION  PTH(lOOl) fGD( 1001) fT' 1001) fC( 1001) fY( 1001) 

DIMENSION  SUBM( 1001 ) fDIAM( 1001  ✓ fSUPM( 1001 ) fBM( 1001 ) 

DIMENSION  SUBT(lOOl)  f DI AT  (1001)  .ScJPTCl 001)  fBT(  1001 ) 

DIMENSION  DIT(1001)fDIM(1001) 

DIMENSION  TEMC1 < 1001 ) »  TEMC2( 1001 >  fTEMC3( 1001 ) 

COMMON/ARE A2/ I MAX  f I HLOW  » I HUP  fDYI fDY2. DY2 , AMB  f  AMT » DT 
COMMON  /AREA1/PA.P0 
DATA  NR/5/ .NW/6/ 

C* 

C*  INPUT  DATA  TO  BE  PROVIDED  FOR  THE  PROGRAM.  IN  THE  READ 

C*  STATEMENTS  'NR'  REFERS  TO  THE  INPUT  DEVICE  NUMBER. 

C* 

READ (NR .99)  MAXSTP.ILIST 

READ ( NR » 99 )  IGRID1 1 IGRID2 

READ ( NR  . 99 )  IGRID3.IPR0F 

READ(NR.IOO)  DELTA»H» DT 

RF AD ( NR *100)  HEATfPAfPO 

READ(NRfIOI)  TINIT fCINIT .PVTOP.PVBOT 

READ (NR .101)  ANUSTTfANUSTB 

READ ( NR v 102)  FIKGt FIKS.FIALFS 

READ (NR. 102)  FIDGf FIDS» ALEW1 

RE AD < NR. 102)  EPS 

99  FORMAT (215) 

100  FORMAT (3F20. 9) 

101  FORMAT  <  4F10.5) 

102  FORMAT ( 3X  .  D9 . 2  f  3X  f  D9  •  2  f  3X  r D9 . 2  > 

C*** 

C***  SAMPLE  DATA  -  (THE  FIRST  AND  LAST  LINES  HERE  INDICATE  THE 

C***  COLUMN  NUMBERS.  THESE  ARE  TO  BE  OMITTED  IN 

C***  THE  ACTUAL  DATA.) 

C*** 

C*** 

C***  123456799012345678901234567890123456789012345678901234567890 


C*** 

c*$* 

4000  100 

C*** 

50  10 

■  1 

c*** 

50  1000 

c**$ 

0.001  0.50 

0.0002 

c*** 

50.000  5.00 

1.0 

c*** 

1.0  0.0  0.5  0.01 

c*** 

200.0  10.0 

c««* 

0.25D+000  5 . 00D+001  0.50D+001 

c*** 

1 . 00D+007  0. 10D+001  1.00D+004 

1 

c*** 

1.00D-001 

c*** 

C***  123456789012345678901234567890123456789012345678901234567890 

C*** 

C*** 

C* 

C*  CREATION  OF  THE  FINITE  DIFFERENCE  GRID 

C* 

CALL  GRIDdGRIDl.  IQRID2fIGRID3fH.D£LTA»Y> 

C* 

C* 
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c** 

c** 

c** 


THE  LOOP  FOR  CALCULTAING  QUANTITIES  AT  NEW  TIME  LEVEL . 

I TIME  =  0 
TIME  =  0.0 
TLOLD  =  TINIT 
TUPOLD*  TINIT 
XLOLD  *  PVBOT 
XUPOLD*  PVTOP 
1 PRINT  =  0 
1000  I TIME  -  I  TIME  +  1 

IF ( ITIME  .GT.  MAXSTP >  GO  TO  1001 
TIME  «  TIME  +  DT 

now  *  T(i) 

TUP  -  T ( I MAX ) 

c* 

C*  NEW  VALUE  OF  THE  SURFACE  COVERAGE  (THETA)  IS  CALCULATED. 

C* 

CALL  CALX ( TLOLD  » TLOW »  XLOLD*  XLOW) 

CALL  CALX ( TUPOLD , TUP » XUPOLD » XUP  > 

CALL  CALTHT  <  TLOW »  XLOW » THETLO) 

CALL  CALTHT (TUP, XUP, THETUP) 

C* 

TLOLD  =  TLOW 
TUPOLD3  TUP 
XLOLD  =  XLOW 
XUPOLD-  XUP 
(PRINT  =  IPRINT  +1 
ICHEK  =  IPRINT/ILIST 

IF  (IPRINT  .NE.  ( ICHEK*  IL  1ST ).  >  GO  TO  555 
WRITE(NW» 149)  TIME 
WRITE ( NW  , 160) 

160  FORMA  T  <  5X  ,  '  TEMP  IN'  ,5X,  'THETA  IN'/) 

WRITE ( NW, 1 50 )  TLOW, THETLO 
WRITE(NW» 161) 

161  FORMAT (5X, 'TEMP  OUT ', 5X ,  ' THETA  OUT'/) 

WRITE(NW, 150)  TUP »  THETUP 

C** 

C**  THE  BULK  CONCENTRATION  OF  VAPOR  IN  THE  SOLID  IS  CALCULATED 

C**  USING  SIMSON'S  RULE-  INTEGRATION. 

C** 

DO  330  I®  1 , IHLOW 

330  TEMCKI)  =  C <  I ) 

DO  331  I = IHLOW, I HUP 

331  TEMC2( I-IGRID1 )  *  C < I > 

DO  332  I  =  IHUP » IMAX 

332  TEMC3< I-IGRID1-IGRID2)  »  C(I) 

11  ~  IHLOW 

12  *  IGRIB2  +  1 

13  *  IHLOW 

CALL  S I MRUL < II * TEMC 1 » D Y 1 » CB I > 

CALL  SI MRUL ( 12,  TEMC2»DY2,CB2) 

CALL  SI MRUL ( 13 » TEMC3 » DY3»CB3 ) 

CBULK  ®  CB1  +  CB2  +  CB3 
WRITE  <  NW# 163 )  CBULK 

163  FORMAT (/I OX , 'CBULK  »  ',F14.8»//) 

C** 

58 


APPENDIX 


c* 

C*  THE  CONCENTRATION  AND  TEMPERATURE  PROFILES  ARE  PRINTED 

C*  ALONG  WITH  THE  VALUE  OF  THE  Y-  COORDINATE. 

C* 

ICHEK  =  I PR I NT /I PROF 

IF < IPRINT  .NE.  <ICHEK*IPROF)  >  GO  TO  555 
URITE(NWf 162) 

162  FORMAT ( 5X » ' CONC « ' » 8X» ' TEMP . ' » 10X  * ' Y  COORD. '/> 

WRITE <NW» 200)  ( C ( I ) f T < I ) » Y ( I ) , 1=1 » IMAX ) 

149  FORMAT <//10X» 'TIME  =  '»F12.8//> 

150  FORMAT  <  3F14 . 6// ) 

C*** 

C***  INCORPORATION  OF  THE  HEAT  TRANSFER  BOUNDARY  CONDITIONS 

C***  AT  THE  OUTSIDE  AND  INSIDE  SURFACES. 

C*** 

555  CALL  TBOUND ( 1 » DY 1 » SUBT » D I AT » SUPT t TLOW * XLOW t ANUSTB »FIKG» DEL ) 
CALL  TBOUND ( IMAX t DY3 » SUBT »DIAT » SUPT , TUP » XUP » ANUSTT » F I K j » DEL ) 
CALL  SOURCE <BT1 »BT2»BT3»FITK»GTB»GTT»T *QTH»BT) 

BT ( 1 )  =  T <  1 )  +  DIAT(l)  -  1  +  SUPT <  1.) 

BT(IMAX)  =  T(IMAX)  +  DIAT < IMAX )  -  1.  +  SUBT(IMAX) 

C*** 

C***  INCORPORATION  OF  THE  SOLUBILITY  CONDITIONS  FOR  VAPOR 

C***  DIFFUSION  INTO  THE  SOLID. 

C*** 

CALL  SOURCE i BM 1 , BM2 » BM3 » FIMNV GMB » GMT  r  C  »  QD . BM ) 

DUM1  =  EPS*THETLO 
BM<1)  =  DUMl/< 1 .+  DUM1) 

DUM2  =  EPS*THETUP 
BM(IMAX)  ■  DUM2/ < 1 ♦ +DUM2 ) 

C#* 

C**  THE  TRI-DIAGONAL  MATRIX  INVERSION  PROCEDURE  FOR  UPDATING 

C**  THE  TEMPERATURE  AND  CONCENTRATION  SOLUTIONS  AT  EACH  TIME 

C**  LEVEL. 

C** 

DO  20  1=1 » IMAX 
DIT  < I )  =  DIAT  < I ) 

20  DIM< I )  =  DIAM<I) 

CALL  TRID< SUBM» DIM»SUPM» BM» IMAX ) 

CALL  TRID<SUBT » DIT » SUPT » BT » IMAX ) 

DO  21  1=1 » IMAX 
C(I)  =  BM(I)  , 

T  < I )  =  BT ( I  ) 

21  CONTINUE 
GO  TO  1000 

200  FORMAT (3F15.S) 

1001  STOP 
END 

C** 

C** 

C** 

C** 
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V 


C<« 

C<« 

C<« 

CX< 

c<« 

c«< 


10 


11 


12 


THIS  SUBROUTINE  COMPUTES  THE  RIGHT-HAND  SIDE  VECTOR  FOR 
TRI-DIAGONAL  MATRIX  INVERSION  PROCEDURE.  THE  TERM 
INVOLVES  TEMPERATURE  (OR  CONCENTRATION)  AT  OLD  TIME 
LEVEL  AND  ALSO  THE  VOLUMETRIC  SOURCE  FOR  HEAT  OR  MASS. 

SUBROUTINE  SOURCE < BM1 , BM2 » BM3 , FIMK »GHB,6MT»C»QD»BM) 

IMPLICIT  REAL<8<A-H,0-Z) 

CQMM0N/AREA2/IMAX, IHLOW. IHUP»DY1 »DY2»DY3»AMB,AMT  ,DT 
DIMENSION  QIK  1001 )  ,BM< 1001 >  ,C< 1001 ) 

DO  10  1=1, IHLOW 

BM( I )  =  C<I)  +  BM1<DY1<DY1<QD< I) 

CONTINUE 

DO  11  I = IHLOW, I HUP 

BM ( I )  =  C(I>  +  BM2<DY2<DY2<QD<  I) 

CONTINUE 

DO  12  I=1HUP , I MAX 

BM ( I  )  =  C(I)  +  BM3<0Y3<DY3<QD<D 

CONTINUE 

BM< IHLOW)  =  C( IHLOW)  +  DT<AMB<AMB<FIMK<QD<IHLOW)/GMB 
BM ( I HUP )  =  C(IHUP)  +  DT<AMT<AMT<FIMK<QD<  IHUP ) /GMT 

RETURN 
END  ' 


C«< 

c<« 


r 

> 

'X 

c> 

C« 

n 

c« 

T 

c« 

‘i 
•  "  •  1 

c« 

% 

A* 

V 

»  • 
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THIS  SUBROUTINE  INCORPORATES  THE  HEAT  TRANSFER  BOUNDARY 
CONDITIONS  ON  THE  OUTSIDE  AND  INSIDE  SURFACES. 

SUBROUTINE  TBQUND< I ,DY, SUBT , DIAT »SUPT,TS,XS» ANUST ,FIKG , DEL ) 
IMPLICIT  REAL<8(A-H,0-Z) 

C0MM0N/AREA2/IMAX, IHLOW, IHUP, DY1 , DY2,DY3, AMB,  AMT,  DT 
COMMON  /AREA1/QA,Q0 

DIMENSION  SUBT< 1001 >,DIAT< 1001 ),SUPT< 1001) 

CALL  THE TAP <  TS  »  XS, DERI V ) 

FI  =  ANUST<FIKG 

F2  =  DEL<DERIV<<  00  +  <QA-QO)<< 1 .-XS)  ) 

DUM1 =  DT/(0.5<DY<DY  -  F2<DY  ) 

DUM2=  <1.+F1<DYXDUM1 


DIAT  < I ) 
SUBT < I ) 
SIJPT(I) 
REfURN 
END 


1.  +  DUM2 

-DUM1 

-DUM1 


C« 

c« 
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(-,#** 

C***  THIS  SUBROUTINE  GENERATES  THE  SPACIAL  GRID  FOR  THE 

C***  FINITE  DIFFERENCE  PROCEDURE. 

C*#* 

SUBROUTINE  GRID< IGRID1 » IGRID2 * IGRID3> H» DELTA* Y > 
IMPLICIT  REAL*8<A-H»0-Z> 

DIMENSION  Y(lOOl) 

COMMON/ ARE A2/ 1 MAX » IHLOW » IHUP ► DY1 , DY2 * DY3 , AMB  *  AMT  ,  DT 

IMAX  ~  IGRID1  +  IGRID2  +  IGRID3  +  1 

IMAX1  =  IMAX  -  1 

I HI GW  =  IGRID1  +1 

IHUP  =  IHLOW  +  IGRID2 

DY1  =  (l.  -  H  -  DELTA ) /DFLOAT ( IGRID1 ) 

DY2  =  DELTA/DFLOAT  < IGRID2) 

DY3  =  H/DFLOAT ( IGRID3 ) 

AMB  =  DY2/IIY1 
AMT  =  DY2/DY3 
Y  < 1 )  =  -1.0 

DO  l  I =2 r IHLOW 

1  Y< I )  -  Y ( I - 1 )  +  DY1 
DO  2  1  =  1 » IGRID2 

2  Y< IHLOW  +  I)  =  Y( IHLOW  +  I  -  1)  +  DY2 
DO  3  1= IHUP , I MAXI 

3  Y  (HI)  =  Y  ( I )  +  DY3 
RETURN 

END 

C*** 

c;*** 


C** 

€**  1  Hi  S  SUBROUTINE  CALCULATES  THE  EXTENT  OF  SURFACE  COVERAGE 

C**  (THETA)  AS  A  FUNCTION  OF  TEMPERATURE  AND  DIMENSIONLESS 

C**  PARTIAL  PRESSURE. 

C** 

SUBROUTINE  CALTHT <  TS » X*  THETA ) 

IMPLICIT  REAL*8< A-H»Q-Z) 

COMMON  /AREAl/QAf QO 
BUM 1  =  <0A  -  QO)/TS 
CONS  =  DEXP(DUMl) 

BUM 2  =  (1.-X)*<1.  -  X  +  CQNS&X ) 

THETA*  C0NS*X/DUM2 

RETURN 

END 

C** 

('■** 


('*** 

f;*x*  THIS  SUBROUTINE  CALCULATES  THE  NEW  DIMENSIONLESS  PARTIAL 

C***  PRESSURE  AT  THE  NEW  TEMPERATURE. 

C**# 

SUBROUTINE  CALX( TSOLDf  TSNEWr  XOLD» XNEW) 

IMPLICIT  REAL*8<A-H»0-Z> 

COMMON  /AREAl/OA.GO 

POW  =  -0 0*<  l./TSOLD  -  l./TSNEW) 

XNEW  *  XGLD#DEXP(POW) 

RETURN 
i  nr 
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r** 

C**  1  HI S  SUBROUTINE  SOLVES  FOR  A  SYSTEH  OF  EQUATIONS  OF  THE 

C,**  FORM  CA3*CX3  =  CB3  r  WHERE  THE  MATRIX  LA3  IS  TR I -DIAGONAL 

c:** 

SUBROUTINE  TRID < SUB f DI AG f SUP f B f N) 

IMPLICIT  REAL*8(A-HfO-Z) 

DIMENSION  SUB< 1001 ) ,DI AG ( 1001 >  f SUP (1001 > ,B(1001 > 

IF (N  .  GT .  1)  GO  TO  10 
B ( 1 )  =  B( 1 )/DIAG( 1 ) 

RETURN 

lv  DO  11  K=2 f N 

RATIO  -SUB<K)/DIAG(K-1 ) 

DIAG(K)  «  DIAG(K)  +  RATI0*SUP(K-1 > 

N  D(K>  -  B < K )  +  RATI0*B(K-1) 

B(N)  =  B(N)/DIAG(N) 

E  N 

DO  1.2  NP1MK=2,N 
K  K  1 

12  B(N  )  -  ( B  (  K  )  -  SUP  (  K )  *B  <  K  + 1 )  ) / D I AG ( K  > 

RETURN 

END 

c** 

i.'.&Xi 

I .  *  -I.  * 

C***  THIS  SUBROUTINE  CALCULATES  THE  COEFFICIENT  MATRIX  [A3 

C***  USED  IN  THE  TRI-DIAGONAL'  FfNITE  DIFFERENCE  PROCEDURE. 

C*** 

SUBROU r I NE  COEFM ( BM 1 f  BM2 f BM3 f F IMK  f  GMB f  GMT  f SUBM ,  D I AM , SUPM ) 
IMPLJ.Cn  REAL*«(A-H»0-Z) 

C0MM0N/AREA2/IMAX f I HLOW  * IHUP  fDYI f DY2 f  DY3 f  AMB f  AMT  ,  DT 
D I MENS 1  ON  SUBM ( 1 00 1 >  f  D I AM ( 100 1 ) f  SUPM (1001) 

CB  =  DT/<DYl*DY1) 

CT  =  DT/(DY3*DY3> 
lMAXt  -  IMAX  -  1 
SUPM  <  1  )  =••■  0.0 
DIAM(l)  -  1.0 
SUBM < 1  )  -  0.0 
SUE M( IMAX)  =0.0 
•D 1  AM  ( I  MAX  )  =  1.0 
SUBM(IMAX)  =  0.0 
DO  1  1=2, IHLOU 
CIJBM(  I)  =  -  BM  1 
SUPM  <  I)  =  -BM1 
J  DIAM(I)  =  1.+  2.*BM1 

DO  2  I "I HLOW f IHUP 
SUBM( I )  =  -BM2 
SUPM< I )  a  -DM2 
DTAM< I  )  =  1,  +  2.*BM2 
DU  3  1 « IHUP f IMAX 1 
SUBM (I)  =  -BM3 
SUPM(I>  =  -BM3 
>  DIAM(I)  =  1.  +  2.*BM3 

SUBM <  ITU.  OW)  =  -2 .  *CB*AMB/GMB 

:»UPM(  IHLOU)  =  ~2.*CB*FIMK/GMB 

DIAM(IHLOW)  =  1.  +  2.*CB*<AMB+FIMK)/GMB 

SUBM (IHUP)  =  -2.*CT*FIMK/GMT 

SUPM< IHUP)  a  -2,*CT*AHT/GMT 

D I AM t IHUP)  *  1.  +  2 . *CT*  <  AMT +F 1 MR ) /GMT 

RETURN 

END 
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c  $  * 

<**  THIS  SUBROUTINE  CALCULATES  THE  DERIVATIVE  OF  THETA  WITH 

C**  RESPECT  TO  TEMPERATURE. 

C** 

SUBROUTINE  THETAP ( TS »  XS »  DERI V  > 

IMPLICI1  REAL *8 ( A— HV 0~Z ) 

XOLH  =  XS 

TSOI.  D-  rs  • 

CALL  CAL.THT  (  TSfiLD *  XOLU  f  THOLD ) 

I SNEW-  t.0001*TS 

BELTS-  TSNEW  -  ISOLD 

CAL  I  CAI.  X  (  TSOLD »  TSNEW » XOLU »  XNEW  > 

CALI  CALTHT (TSNEW. XNEW ,THNEW) 

DERIV  I  JHNEW  -  THOLD ) /BELTS 

RE i URN 

END 

!  Ek 
I  .  .*  K  ' 


cfcX*  THIS  SUBROUTINE  CALCULATES  THE  INTEGRAL  OF  A  FUNCTION  USING 

Ck*:*  SIMOON'S  RULE 

(.;*** 

SUBROUTINE  SIMRUL I N p G » DX t SIMS ) 

IMPLICIT  REALMS ( A- H » Q-Z ) 

DIMENSION  GtlOOl) 

NK  N-J. 

NN  N -2. 

01  -  G < 1 )  +  GIN) 

c»;>  -  o.o 

DO  to  K~2 » NK » 2 
Jo  02  -  G2+  OIK) 

0.5  :;=  0.0  .  • 

JPINN  .EQ.  1)  GO  TO  15 
DO  1!  K~3»NN»2 
.1.1  Ho  -■  <t>3  +  GIK) 

IS  SIMS  ='  DX*IGl+4.0*G2+2.0*G3>/3.0 

RETURN  • 

END 

END  Of  FILE' 
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